#ifndef ADAMS_MOULTON_HPP
#define ADAMS_MOULTON_HPP

#include "AdamsMoultonCoeff.hpp"
#include "LMM.hpp"
#include <cmath>
#include <iostream>

// constexpr int ITE = 20;

class AdamsMoulton : public LMM {
public:
  AdamsMoulton(int dim, double dt, double t0, double t1, int p)
      : LMM(dim, dt, t0, t1, p) {
    if (p > IntScheme::AMCoeff::maxP || p < IntScheme::AMCoeff::minP) {
      std::cerr << "Invalid p" << std::endl;
    }
    m_s = IntScheme::AMCoeff::s[p - IntScheme::AMCoeff::minP];
    m_u.resize(std::floor((m_t1 - m_t0) / dt) + 1, m_dim);
    m_u.setZero();
  }

  const std::vector<double> getAlpha() const {
    int idx = IntScheme::AMCoeff::offset[m_p - IntScheme::AMCoeff::minP];
    return std::vector<double>(IntScheme::AMCoeff::alpha.begin() + idx,
                               IntScheme::AMCoeff::alpha.begin() + idx + m_s +
                                   1);
  }
  const std::vector<double> getBeta() const {
    int idx = IntScheme::AMCoeff::offset[m_p - IntScheme::AMCoeff::minP];
    return std::vector<double>(IntScheme::AMCoeff::beta.begin() + idx,
                               IntScheme::AMCoeff::beta.begin() + idx + m_s +
                                   1);
  }

  int getS() const { return m_s; }
  Eigen::MatrixXd getU() const { return m_u; }
};

#endif // ADAMS_MOULTON_HPP